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Abstract 



We consider the recovery of real- valued bandlimited functions from the absolute values of 

i— i their samples, possibly spaced nonuniformly. We show that such a reconstruction is always 

<1^ possible if the function is sampled at more than twice its Nyquist rate, and may not necessarily 

be possible if the samples are taken at less than twice the Nyquist rate. In the case of uniform 

f-H samples, we also describe an FFT-based algorithm to perform the reconstruction. We prove that 

it converges exponentially rapidly in the number of samples used and examine its numerical 

r* behavior on some test cases. 
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o 1 Introduction 

o 

In a series of recent papers 03 El El, Balan, Casazza and Edidin have investigated the possibility of 
reconstructing finite-dimensional signals using measurements that do not contain sign or phase in- 
formation. Motivated by an application in the denoising of speech signals, they studied M-element 
frames for W 1 , i.e. collections of M vectors that span W 1 . They considered M-element frames, 
{fk}i<k<M> fk £ ^"» sucn that any vector xeK" can be uniquely determined from the inner prod- 
ucts {\(fk,x)\}i<k<M U P t° an ambiguity of a sign factor. Using frame theory and combinatorial 
methods, they showed that such frames exist if and only if M > 2n + 1. In [1J, a computational 

n(ji I ]\ 

method to carry out this reconstruction was described in the case where M > v 2 - , using a special 
class of frames. 

It is natural to ask if there are analogous results for continuous-domain signals, namely in 
the context of samples of bandlimited functions. The well-known Whittaker-Shannon-Kotelnikov 
(WSK) sampling theorem [4] shows that if a bandlimited function / is sampled at a rate greater 
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than its Nyquist frequency, then it can be uniquely reconstructed from the samples. However, if 
the signs of the samples are unknown, this condition may no longer be sufficient. An example 
illustrating this is given by the functions sin (n [z + \) ) and cos (n (z + 1)), which agree in abso- 
lute value at z = §, k £ Z. On the other hand, it might be expected that if we oversample / at a 
sufficiently high rate, the samples may contain enough redundancy that can we afford to lose their 
signs and still recover / up to a global sign factor. It turns out that this is indeed the case. 

In this paper, we use a complex variable approach to show that if a real-valued bandlimited 
function / is sampled at more than twice its Nyquist rate, then / can be uniquely determined from 
the absolute values of its samples up to a sign factor. Conversely, we find that if / is sampled at 
less than twice its Nyquist rate, then it is not always possible to uniquely determine it in this way. 
We present an algorithm to perform this reconstruction, and show that it converges exponentially 
rapidly in the number of samples used. We consider a fairly general class of nonuniformly spaced 
samples in this paper, although our numerical approach is developed with uniformly spaced sam- 
ples in mind for reasons of computational efficiency. 

We review some existing theory on nonuniform sampling and bandlimited functions in Section 
[2j and then state and prove our main theoretical results in Section |3j We describe our algorithm 
and study its convergence properties in Section |4j and apply it to two test cases in Section [5} 

2 Background Material 

We normalize the Fourier transform as /(ft)) = §™^f{t)e~ l7Zim dt for Schwartz functions / and 
extend it to tempered distributions in the usual way. For < p < °°, we define the Paley-Wiener 
spaces of bandlimited functions [14J by 

PW/ = {/ 6 Z/:supp(/)c [-!,!]}. 

An entire function g is said to be of exponential type b if 

fc = inf(jB:|s(z)|<ePM,zeC), 

and we denote this by writing type(g) = b. By the Paley-Wiener-Schwartz theorem |[8), PW£ can 
be equivalently described as the space of all entire functions / with type(/) < Kb whose restric- 
tions to R are in LP . It also follows that PWi C PW? for p < q. Functions / G PWE satisfy the 
classical estimates H/'H^p < Tib \\f\\jj> and ||/(- + ic)\\ LP < e Kb ^ ||/|| L p, respectively known as the 
Bernstein and Plancherel-Polya inequalities lfTTl[T4ll . 

We now consider a sequence of points X = {xt] C K, indexed so that Xk < Xk+\- For any set 
B, we denote the number of Xk in B by N(X , B) . We say that X is separated if infyt \xk+ 1 — JCjt | > 0. 
Following [Q, X is also said to be uniformly dense if it satisfies 



sup 

k 



k 

Xk-~; 
d 



<oo (1) 



for some finite d > 0. We denote this by writing D(X) = d, and by D(X) = °° if X is not uniformly 
dense. We will mainly deal with separated, uniformly dense sequences in this paper. It is worth 



mentioning that D(X) is not directly related to Beurling's upper and lower densities (see IIT410 . and 
there are sequences X with finite Beurling densities but for which D(X) = °°. 

The generating function of a sequence X C R is given by 

S(z)=z 5x hm J] i 1 --), (2) 

'^o<N<A x k J 

where 8x = 1 if G X and 5x = otherwise. For a uniform sequence^ = |, 5(z) = sm ^ z ^ . If X is 
separated and D(X) = b <°°, then the limit in ([2]) is finite and the function S lies in the Cartwright 
class CWi, O, the set of all entire functions / with type(/) < Kb that satisfy the growth condition 

max(log |/(Q |,0) ^ 
» — at < oo. 

t 2 + l 

In fact, functions in CW\, satisfy the apparently stronger condition |[T0l 

|log|/(* 



t 2 +l 



dt < oo. (3) 



The following result shows that an arbitrary function in CWt, can be expanded in the form ([2]) and 
gives a useful geometric description of its zeros. [fTTTl 

Theorem. (Cartwright-Levinson) Let W + (6 ,r), W~(0,r) andW'(6,r) respectively be the wedges 
{z: \z\ <r, | argz| < 0}, {z: \z\ <r, |tt — argz| < 6}and{z: \z\ < r, |argz| > 0,\n-argz\ > 9}. Sup- 
pose f E CWt,, f ^ 0, and let U = {u^} be the set of its zeros. Define b' — ^ inf(type(e ,COz /(z))> CO G 



1: For any 6 G (0, |), as r — > oo, 

N(U,W + (6,r)) b' N(U,W-(6,r)) b' N(U,W'(6,r)) 

> 7T; > Wi an d > 0. 

r 2 r 2 r 



2: U satisfies £i M , i >0 IMi) I < °° 



«jtl><J I y u k - 



for some constant q. 



f(z)=f(0)e^z 8u \im J] fl--V (4) 

r ^~o<L|<A u k J 



of / satisfies lim^oo - r ' < £. The expression on the right side of M is called a canonical 



If / is real-valued on K, then b' = b and (7 = 0, so in particular, the sequence of real zeros Vet/ 
of / satis 

We will also need a second, deeper theorem on CWt,. 



Theorem. (Beurling-Malliavin) Let f G CWy. Then for any e > 0, there exists h G CW E , h ^ 0, 
such that fh G L°°. 

A discussion of this result and its significance can be found in [7J. We can actually choose the 
above h so that fh G LP for any p > 0, by replacing h with h(j) where is a function for which <p 
is smooth and has sufficiently small support. By considering the function h(z)h(z), h can also be 
taken to be real and nonnegative. 

For 1 < p < °°, X is called a sampling sequence for PWF if there are constants a p ^ and j5 p j, 
such that for all / G PW^, a Pib \\f{X)\\ lP < \\f\\ LP < Pp,b\\f( x )\\iP- In particular, this condition 
implies that any / G PWF is uniquely determined by its samples at X. Precise geometric charac- 
terizations of sampling sequences can be very complicated (see l|4l [T3l[T4| ). but for real, separated 
and uniformly dense sequences X, it is necessary that D(X) > b and sufficient that D(X) > b for 
X to be a sampling sequence. If D(X) > b and S is the generating function of X, a consequence 
of the Beurling-Malliavin theorem is that there is an h G CW £ with zero set U, £ = D(X) — b, such 
that any / G PW£° can be expressed in terms of its samples, 

/(!) % E J W WWM)' <5) 

with uniform convergence on compact subsets of C\U [|9l [121. -D(X) can be thought of as a 
generalization of the "sampling rate" to uniformly dense sequences X, and when X = {-} is a 
uniform sequence, the condition b < D(X) = s simply says that / is being oversampled beyond its 
Nyquist rate, h can be taken as constant for uniform sequences, and the expansion ([5]) reduces to 
the classical WSK sampling theorem. 

3 Main Results 

We fix p G (0, oo] for the rest of this section. We first show that if a bandlimited function is sampled 
at more than twice its Nyquist rate, then we can reconstruct it up to a sign factor from the absolute 
values of the samples. 

Theorem 1. Let f G PW^ be real-valued on M, and let X C K be a separated, uniformly dense 
sequence with D(X) > 2b. Then f can be uniquely determined from a& = \f(xk) \, up to a sign factor. 

Proof. We normalize b = 1 without loss of generality. Let S be the generating function of X and 
suppose h is given as in ((51). The zeros of / and h are countable, so we can choose c > so that / 

— /I 

and h have no zeros on the line L = {z : Im(z) = c}. We can write g — /(• + ic) 2 G PW P as 

n = y 2 S(z + ic)h(x k ) 

[Z) k hS k Hz + ic)S'(x k ) (z + ic - x k ) ■ W 

From Bernstein's inequality, g' G PWJ? and differentiating ( 6 ) gives a similar expansion for g'(z) . 
Let a* k = f(x k + ic) and choose a point xi G X. Since / has no zeros on L, there is a branch of 
arg/, which we denote by arg /, that is continuous on L and satisfies arg /(jc/ + ic) G (— ic,7c]. 
We define arg ^ in the same way, and we then have 



a k = ex p( 2 lo § \s( x k)\ + iarg f(x k + ic) j (7) 

= ^^(^l^exp^^Im^^r + argo^)^ (8) 

where T] = exp (i Qarg g(;c/) — arg /(jc/ + /c))) = ±1. We have now determined samples of /, 
which we can use to express / as 

f( z )=n f at 5(z ~ /c)/ifa) (9) 

D 

Remark. The proof of Theorem [T] suggests a three- step procedure to recover / from {a k }. We can 
determine f 2 from {a k } and take its square root by unwrapping its phase. Since / will typically 
have zeros on the real axis, we first move up in the complex plane using Q, unwrap the phase 
there with ([8]) and then move back to the real axis with Q. We will use this approach in Section [4] 

The next result shows that Theorem [T] is in a sense sharp. If we sample a function at less than 
twice its Nyquist rate, it may or may not be uniquely determined by the absolute values of the 
samples, essentially depending on "how many" of the samples are zero. 

Theorem 2. Let X C M. be a separated sequence with D(X) < lb. Then there is a real-valued 
function f E PWE that cannot be uniquely determined from \f(x k ) \upto a sign factor. If in addition 
D(X) > %b, then there is another real-valued function f E PW£ that can be uniquely determined 
from \f(xk)\ up to a sign factor, but for which there is no subsequence Y C X with f(y k ) = and 



Proof As before, we normalize b — 1. Suppose D(X) = 2 — e for some £ > 0, and let X\ = {x2 k } 
and X 2 = {x 2k+ i}. It follows from the definition (TH) that D(Xi) = D(X 2 ) = 1 - f . Now let 5i 
and S2 be the generating functions of X\ and X 2 . By the Beurling-Malliavin theorem, we can find 
h\,h 2 G CW £ / 2 such that S\h\ and 52^2 are in PWf. Then the functions f\ = S\h\ + S 2 h 2 and 
f 2 =Sih -S 2 h 2 are also in PWf and satisfy \fi(X)\ = \f 2 (X)\. 

For the other part of Theorem 2 suppose £ above satisfies £ < |. By countability, we can 
choose c so that h-^iz) = h\{z + c) has no zeros on X. Let / = Si/13, so that among the samples 
of / at X, only the ones at X\ are zero. If there is any real- valued g E PWf with \f(X) \ = \g(X)\, 

then f 2 — g 2 is in PW 2 P and has zeros at X, so by considering canonical product expansions, 
f 2 — g 2 = S\S 2 h4 for some h\ ECW £ . Since g 2 = S[(S{h 2 — S 2 ha) has an analytic square root, it can 
only have double zeros, so in particular, S\h 2 — S 2 h^ has to be zero on X\ . This implies that h.4 must 
have zeros onX\, which contradicts the Cartwright-Levinson theorem because D(X\) — 1 — f > £. 
So h 4 = and f = ±g. □ 

Remark. The nonexistence of Y in Theorem[2]is what makes the result interesting, as it means that 
/ cannot be determined from Y by just using a canonical product. In other words, the nonzero 
samples at X 2 play a role in the uniqueness of /. However, there appears to be no simple charac- 
terization of all such functions / or a numerically useful method of computing / from \f(X) \ . 



Remark. We have not considered the border case of D(X) = 2b in the above results, in which case 
the conditions required on the sequence X would become more subtle and depend on the value of 
p. However, in the elementary case where Xk = A is a uniform sequence and p = 2, the conclusion 
of Theorem [T] still holds by just using the WSK sampling theorem in place of ([6]) and Q. 

There are no simple analogs of these results if we allow / G PW£ to be complex- valued. In 
Theorem [I] we used the fact that when / is real- valued, f 2 has the same samples as \f\ 2 , but this 
is no longer the case for complex-valued /. In general, such an / will have complex zeros Uk and 
complex- valued functions of the form Bf, where B is a Blaschke product formed from any subset 
of {Uk}, will be in FWF and satisfy \Bf\ = \f\ identically on IR. If we require all the zeros of / 
to be real, then since / can be written as a canonical product over them, it is simply a modulation 
of a real- valued function g, i.e. it has the form f(z) = e'^ cz+d 'g(z) for c,(/£l This situation 
is in contrast to the findings in ||2), where the types of frames the authors studied had results for 
complex vectors comparable to those outlined in Section [T]for real vectors. 

4 A Reconstruction Algorithm 

We now describe how to computationally implement the technique in the proof of Theorem [T] 
We restrict our attention to uniform sampling sequences here, as they lead to convolution-type 
sampling series that can be calculated efficiently by Fast Fourier Transform (FFT) methods, but 
the same ideas can be adapted to the nonuniform case. We first define the following functions: 



G(z,M) = s^le- 

nz 



I'M- 



cos(7Tz) sin(7Tz) sin(7Tz)\ ?l z 2 



/, . / WJI/K.I MiM'^l ^"I'KI \ 

G(zM) = — - — - — -)e 

y ! \ nz M nz ! 

1 f z 
G*{z,M) = — / G(t,M)dt 



Mj z - M 

We also denote the strip {z : |Im(z)| < 8} by T§. Then we have the following results from [12J, 
reproduced here in slightly different forms for our purposes. 

Theorem 3. (Schmeisser-Stenger, 2002) Let f G PW^, s > b and d < 1. Then for Re(z) G 
[-dM,dM], f^)=Lt-Mf( K s )G(z-k 1 M)+E(z), where \E(z)\ <CiM-VV^(H)^+2*M Z )l , 
L oo for some constant C\ = C\ ( -). 



Theorem 4. (Schmeisser-Stenger, 2002) Suppose f is analytic in T§ and \f(z)\ < K\z\ for zET§ 
and some constant K. Then for z G R,-| < Z < f , f(z) = li^ dM f{ji) G i Mz ~ k ^ M ) + E ( z )> 
where \E(z)\ <C2K(M/8) l / 2 e ~ \\f\\ L °° for some constant C2. 



Theorem [3] is a form of (|5]) with a non-bandlimited h. It converges very rapidly in practice for 
about s > 1. 3b. Theorem [4] is a version of Theorem [3] as s — V °°. 

These results lead to the following numerical approach. We impose a mild restriction to rule 
out some pathological functions, and assume for notational convenience that an odd number of 
samples are used. 

Algorithm 5. Suppose f G PWjf is nonzero in some strip {z : |Im(z) — c| < 8}, 8 > 0. Let s > 2b 
anda k = \f(*)\for -M <k<M. 

1: Compute g M (§) = I%L_ M a 2 k G(z-k+ic,M) atz= ft, -M 2 <n<M 2 . 

2: Compute g' M (^) =Lf=-M a l G '( z ~k + ic,M) for z as above. 

3: Compute Q M (n) = il^^LgJm [ ^4?V I G*(Mn-k,M) for -(M-2) < n < M- 1. 

\gM{ Ws )J 

4: Compute /?m(«) given by 7?m(0) = 0, #m(«) = ^m(« - 1) + <2m(" - 1) for -(M - 1) < n < 
M-\. 

5: Compute / M (f)=lf^ (M _ 1) V / |^(5)l^ ( ^ W+argSM(0)) G(z-fc-/c,M). 

The convolutions in steps 1 and 2 can be performed by 2M FFTs of size 2M + 1 each. Steps 3 
and 4 are most efficiently done by direct computations that respectively involve 0(M 2 ) and O(M) 
operations. Step 5 involves 2N FFTs of size 2M+ 1 for some integer N, depending on how finely 
we want to compute /. This gives an overall complexity of 0(M 2 logM). 

Theorem[4]is used above to calculate the integral in ([8]). The advantage of this approach is that 
the error bound in Theorem [4] does not depend on the derivatives of /. If we instead used con- 
ventional methods of integration such as Gauss quadrature, we would need additional restrictions 
on 8 to ensure that the algorithm actually converges. There is no simple closed-form expression 
for G* but its numerical calculation only depends on M, so we can tabulate its values by a stan- 
dard quadrature method using 0(M 2 ) operations and reuse them for different / with a lookup table. 

We can establish the following convergence result for Algorithm [5J 

Theorem 6. Let f and /m be defined as in Algorithm [5] As M -> °°, /m -^ Vf uniformly on 
compact subsets of R, where 7] = ±1 and \f M (z) - r)f(z)\ = O ^-"KkC 1 "? ),f)M+4xc\ 

This estimate is somewhat conservative and larger exponents of convergence are possible if 
we make more assumptions on /, but it shows how 8 and c affect the rate of convergence. Many 
real- world bandlimited signals have most of their zeros on or near the real axis, so choosing c too 
small will result in a small 8 while choosing c too large will sharply increase the constant in the 
error bound. Values of c between about 0.01 to 0.25 appear to work well in practice. 

Before we prove Theorem [6l we will need an auxiliary lemma. 

Lemma 7. Let h e PW^ be nonzero in T§. Then -^l < K{ |z| 2 + 1 ) in T§/ 2 f or some constant K. 



Proof. Let U = {u k } be the zeros of h(z — id) lying in the upper half plane C + = {z : Im(z) > 0}, 
and suppose £ < Im(z) < 4£ . The Plancherel-Polya inequality shows that h{- — id) G PW|° and that 
the function e mbz h(z — id) is bounded and analytic on C + , so it has the inner-outer factorization 

m 

logh(z-i8) = ir \+f ]og\h(t-iS)\dt + ]og(ll lu l^ llz - U l \+Az + B 



for some constants A and B. We differentiate this to find that 



z-u k 



h'(z-id) 



h(z-iS) 



1 
Ki 



log |/i(? — id) | ^ 2flm(wfc) 
"2 at + 2^ 



it-zf 



(z-u k )(z-u k ) 



+A 



< (2/5) 2 (|Re(z)| 2 + l) 



lo g \h(t-i8)\ dt+ ^rn(u k ) +lAll 



t 2 + l 



kl 2 +i 



Since PW^ C CWJ,, condition ([3]) shows that the integral in the first term is finite. The Cartwright- 
Levinson theorem implies that the sum in the second term also converges, which finishes the proof. 

□ 

It is possible to obtain sharper results than this, but this lemma is good enough for our purposes. 



Proof of Theorem^ We proceed by establishing several intermediate bounds and then combine 
them all at the end. To simplify the notation, we will always use z EM. and k G Z to denote func- 
tion arguments within norms, e.g. the norm of F G L°°(IR) will be denoted by ||F(z)||. 



Let g and r\ be defined as in the proof of Theorem [Tl Define I r M = [ — ( [r(M — 1)J — 1) , [r(M 
1)J — 1] . i r ,u = %P\Ir,M and j rM = %f) [-M[rM\,M[rM\]. Theorem g shows that 



\sM(m- 



< 



Ci 



e -f(l-¥)M+2*c ||f|| 2 



The error term E in Theorem[3]is an entire function and it satisfies the classical Cauchy estimate 

\\E'(z + ic) 



lz ( , 4u :1 <27rmaxf \E(z + ic)\, tz: |Re(z)| < \(M+ l),jlm(z)| < -L 



This shows that 



^2 :— ||#m(m^ 



,// k 



-Ms!\\l°°{h /iM ) 



< 



7 ^l e -% (1-f )(M-l)+2*c+l || f ||2 



M 



A standard identity for Jacobi theta functions (" IU5II . p. 475) gives the bound 
M 



E 3 :- 



£ G{z-l-ic,M) 

l=-M 



< 



L°°{h/ija) 



£ e &(«M*-o 2 ) 



L°°(R) 



< - e iM V2M, 

~~ 2 



and similarly, 



E 4 :-- 



(k+l)M 

£ \G*(Mk-l,M)\ 

l=(k-2)M 



l °°ih/4 



< -V2M. 
~ 2 



M+2) 



Now by Lemma 7 |gi < K(\z\ 2 + 1) in T S / 2 for a constant K. Since \\g M -^IIl~(/ 3/4m+2 ) "^ ° 



and g has no zeros on E, we have || 1/gM ||/w/ ) < °° for sufficiently large M. We also have 

£ bound 

8m\Ms) 8 \Ws> 



8m(w s ) s(m) 



l°°(h/A,. 



M I 



< 



< 



1 



8M 

1 



L"(l 



E 2+\W\\ L ~(I. 



3/4.M+2! 



ih/AM+l) 



8m 



L-ih/4. 



M+21 



gM 



L°°(l 



(K(M 2 + l) + l)max(E h E 2 ). 



3/4.M+2J 



We can use Theorem |4] with this to find that 



R M (k) 



Mm 



clt 



< M 



< M 



< 



Q M (k)- / Iml^fe \dt 

k-\ 



8m\m s ) 8 \Ms 



gM{ Ws ) g( m ) 



(*3/4,m) 

^fe/4,M\{l-Ll(M-l)J}) 

£ 4 + C 2 (2M/5) 1/ V^ 



° (h/i t 



M I 



> \Ms< 



'\Ms) 



gM 



L°°(h/4M+2) 



(K(M 2 + 1) + 1) M 3/2 E 2 + C 2 KM 3/2 (2/8y 1/2 e- 



'U3/4,m), 
k8M 



We now put everything together. Using the elementary inequality | \u\ l > 2 e'° — |v| l ' 2 \ < \u — v\ l ' 2 + 
1 1 — e ld | |v| l ' 2 along with M gives 



/«G)-n/ 



< 



^°°(A/4,m) 



£ 3+ C ie -f( 1 -?)^- 1 ) +27rc 



/-a 



3/4,m; 



+ 



V e -2(^T)(z-') 2 +^' 



< 



k 



L°°(R) 



/m ( - + «c 



gM 



1/2 /, 



1°°('\,m\'3/4,m) 

E 3 



"('3/4,1 



+2Cie-t(i-f)(^-D+4^ 



L°°(R) 



< \e\ /2 + 



1 -exp U M (k) - f Im (^rPj ^ + argg M (0)-argg(0) 

^ W £3)+2C 1 e-f( 1 -f)^ 1 H 4 -|| / || L . (R) 

= C3 (/ )M 7 / 2 e X p (- mi n ( * (l - f ) , f ) M + max g,4,c 
for some C3(/) < °° and sufficiently large M, which establishes the result. 



'°°('3/4,m) 



□ 



5 Numerical Experiments 

We illustrate how Algorithm [5] works on two test cases. We consider the translated Bessel function 
(see [15J) f(z) =Ji(z + 20) sampled at z = k, —M <k<M, and a collection of 2M + 1 samples 
taken from an 8-bit, 44khz audio file. To show the effect of the parameter c, we take c = 0. 1 in 
the first example and c = 0.04 in the second one. Values of G* are tabulated using the built-in 
Gauss-Kronrod algorithm in MATLAB. We measure the worst-case error over I\/2,m+\ as defined 
in Section HI in order to avoid influence from inaccuracies close to the boundary of the full domain 
h,M+\- 



M 


Error over /j^^+i 
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3.7490- 10- 2 


20 


5.9513- 10 4 


30 


4.0158 10" 5 


40 


3.8732- 10" 6 


50 


3.8362 -10- 7 




-50 -40 -30 -20 



Figure 1 : The graph of f(z) = J{(z + 20) and the reconstruction errors with c = 0. 1 . It is clear that 
the error decays exponentially as M increases. 
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Error over 7 1/2; m+i 
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5.1622- 10" 3 
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40 


1.1563- 10 4 


50 


8.8637-10 " 5 




Figure 2: The graph of a section of audio data and the resulting reconstruction errors with c = 0.04. 
The convergence here is slower than it was in the preceding example. This is likely due to the 
presence of complex zeros with imaginary parts very close to c, as well as the proximity of the real 
zeros. 
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